module normeanvar
!Mean and Variance of the moment conditions of sup-test 
!based on the complete vector of v_eta
implicit none
contains
subroutine meanvar(matA,matB,MD,VD,b,eta,si,sigma,dmut,dvsig,n,p)
use testmoments
use testmoments2
use linear_operators
use ghs
use matrixfuns
use datahadler
integer, intent(in)::n,p
double precision, intent(in)::b(n),eta,si,dmut(p,n),dvsig(p,n*n),sigma(n,n)
double precision, intent(out)::matA(p,p),matB(n+1,p),MD(p+n+1),VD(p+n+1,p+n+1)
double precision:: isigma(n,n),rsigma(n,n),irsigma(n,n),beta(n),vejmem(1,n),dmach,tol
integer i,j,irank

isigma=.i.sigma
tol = 100.0*dmach(4)
call dchfac(n,sigma,n,tol,irank,rsigma,n)
rsigma=.t.rsigma
!rsigma=.t.chol(sigma)
irsigma=.i.rsigma
beta=matmul(.t.rsigma,b)
forall(i=1:p)
MD(i)=0.0d0
end forall
MD(p+1)=.25d0*ejm2(beta,eta,si)-.25d0*(n+2.0d0)*n
MD(p+2:p+n+1)=matmul(rsigma,ejmem(beta,eta,si))

matA=matmul(dmut,matmul(isigma,.t.dmut))&
&+.5d0*matmul(dvsig,matmul(kron(isigma,isigma),.t.dvsig))



VD(1:p,1:p)=matmul(dmut,matmul(isigma,.t.dmut))&
&+.25d0*(((dvsig.x.kron(.t.irsigma,.t.irsigma))&
&.x.(emempkememp(beta,eta,si)-(vec(dble(eye(n))).xt.vec(dble(eye(n))))))&
&.x.(kron(irsigma,irsigma).xt.dvsig))&
&+.5d0*((dvsig.xt.kron(irsigma,irsigma)).x.(emkememp(beta,eta,si).x.(irsigma.xt.dmut)))&
&+.5d0*(.t.(((dvsig.xt.kron(irsigma,irsigma)).x.(emkememp(beta,eta,si).x.(irsigma.xt.dmut)))))

VD(p+1,p+1)=(ejm4(beta,eta,si)/16.0d0)+((n+2.0d0)*(n+2.0d0)/4.0d0)*ejm2(beta,eta,si) &
&+(n*n*(n+2.0d0)*(n+2.0d0)/16.0d0) -((n+2.0d0)/4.0d0)*ejm3(beta,eta,si)&
&+(n*(n+2.0d0)/8.0d0)*ejm2(beta,eta,si)-(n*n*(n+2)*(n+2)/4.0d0)-MD(p+1)*MD(p+1)


VD(p+2:p+n+1,p+2:p+n+1)=ejm2ememp(beta,eta,si)+(n+2.0d0)*(n+2.0d0)*eye(n)&
&-2.0d0*(n+2.0d0)*ejmememp(beta,eta,si)
VD(p+2:p+n+1,p+2:p+n+1)=(rsigma.x.VD(p+2:p+n+1,p+2:p+n+1)).xt.rsigma

forall(i=p+2:p+n+1,j=p+2:p+n+1)
VD(i,j)=VD(i,j)-MD(i)*MD(j)
end forall

VD(1:p,p+1)=matmul(dmut.xt.irsigma,.25d0*ejm2em(beta,eta,si)-.5d0*(n+2.0d0)*ejmem(beta,eta,si))&
&+.5d0*matmul(dvsig.xt.kron(irsigma,irsigma),vec2(.25d0*ejm2ememp(beta,eta,si)&
&-.5d0*(n+2.0d0)*ejmememp(beta,eta,si)-eye(n)*(.25d0*ejm2(beta,eta,si)-.5d0*n*(n+2.0d0))))

vejmem(1,:)=ejmem(beta,eta,si)
VD(1:p,p+2:p+n+1)=(((dmut.xt.irsigma).x.(ejmememp(beta,eta,si)-(n+2)*eye(n)))&
&+.5d0*((dvsig.xt.kron(irsigma,irsigma)).x.(jimemkememp(beta,eta,si)&
&-(n+2.0d0)*emkememp(beta,eta,si)-(vec(dble(eye(n))).x.vejmem)))).xt.rsigma

VD(p+1,p+2:p+n+1)=matmul(rsigma,.25d0*ejm3em(beta,eta,si)-.75d0*(n+2.0d0)*ejm2em(beta,eta,si)&
&+.25d0*(n+2.0d0)*(3.0d0*n+4.0d0)*ejmem(beta,eta,si))

forall(i=p+1:p+n+1,j=1:p)
VD(i,j)=VD(j,i)
end forall
forall(i=p+2:p+n+1)
VD(i,p+1)=VD(p+1,i)
end forall
matB(1,:)=-matmul(dmut.xt.irsigma,ejmem(beta,eta,si))&
&-matmul(dvsig.xt.kron(irsigma,irsigma),.5d0*vec2(ejmememp(beta,eta,si))&
&-.5d0*(n+2.0d0)*vec2(dble(eye(n))))

matB(2:n+1,:)=-(rsigma.xt.emkememp(beta,eta,si)).x.(kron(irsigma,irsigma).xt.dvsig)

end subroutine meanvar
end module normeanvar